Identification of optimal contraction sequences for tensor networks 
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Tensor network Ansatze provide powerful tools for the study of quantum many-body systems on a 
lattice in the low-energy regime, representing the state of a system as an efhciently-contractible net- 
work of multi-index tensors optimised numerically by means of a variational algorithm. The efflcient 
contraction of tensor networks is vital to both the development and the implementation of tensor 
network algorithms, but determination of optimal contraction sequences is presently performed by 
human operators in a process which is both tedious and prone to error. This paper presents an 
algorithm for exhaustively searching the space of contraction sequences on a practical timescale, 
automating the determination of both the optimal contraction sequence and the cost function as- 
sociated with the contraction of a given tensor network. By removing the burden of computing 
these contraction sequences from the researcher, this tool both facilitates the development of novel 
efhciently-contractible tensor network Ansatze and assists in the optimal implementation of existing 
tensor network algorithms. 

PACS numbers: 02.70.-c, 05.30.-d, 03.67.-a 
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I. INTRODUCTION 

Tensor network algorithms provide powerful tools for 
the study of a wide variety of physical systems. They are 
perhaps best known for their use in condensed matter 
physics as numerical techniques for the study of quantum 
many-body systems on a lattice [e.g. [l|-|24|. but recent 
breakthroughs blending ideas from quantum information 
with advanced numerical techniques have led to the con- 
struction of new Ansatze [e.g. l25l - [28j having applications 
in fields as diverse as holography and the AdS/CFT cor- 
respondence |29l - [3l| and the classification of topological 
phases in quantum spin systems [32h35| |. As a numeri- 
cal tool, a tensor network algorithm typically comprises 
an Ansatz for the description of pure or mixed quantum 
states, which is composed of a network of tensors, and 
an iterative procedure for updating this Ansatz. Exam- 
ples include the Density Matrix Renormalisation Group 
(DMRG) m, 113 and Time Evolving Block Decimation 
(TEBD) algorithms [Mlli, both of which are based on 
the Matrix Product State (MPS) Ansatz, and also Tree 
Tensor Networks (TTNs) (t^l, Projected Entangled Pair 
States (PEPS) [M Ell, and the Muhi-scale Entangle- 
ment Renormahsation Ansatz (MERA) [25l - l28| . 

The fundamental challenge to the numerical study of 
quantum many-body systems on a lattice is that the 
number of degrees of freedom, and thus the computa- 
tional cost associated with exact simulation, grows expo- 
nentially with the size of the system. To overcome this 
challenge, tensor network Ansatze replace the coefficients 
Cii...i„ of a quantum state 1-0) 
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with a network of tensors whose dimensions are such that 
the number of coefficients required to describe the net- 
work exhibits a better scaling in n, the number of lat- 
tice sites, than does the number of coefficients Ci^,,,i„ in 
Eq. ([T]). Indeed, for many tensor networks this scaling in 
n is polynomial rather than exponential. 

Given this reduction in the number of coefficients in 
the description, a tensor network Ansatz is capable only 
of representing states which lie within some restricted 
region of the Hilbert space of the system, but neverthe- 
less these Ansatze and associated algorithms are capa- 
ble of providing substantial insight into the physics of a 
wide variety of systems in appropriate regimes (again see 
Refs. [l|-|24l for examples). 

In order for a tensor network algorithm to be useful 
as a tool for numerical computation, it must be possible 
to perform the operations of the algorithm for a reason- 
able computational cost. An economical description of 
the relevant part of the Hilbert space of a system is a 
good start, but this is not the only factor which must be 
taken into account: When determining whether a given 
tensor network algorithm is computationally feasible, the 
structure of the tensor network itself also plays a signif- 
icant role. In describing scaling of the cost of a tensor 
network algorithm, it is customary to express this in the 
form of a polynomial in some refinement parameter Xy 
which may (for example) correspond to the dimensions 
of indices within the tensor network. Assuming that the 
coefficients of this cost polynomial are small, costs which 
scale as excessively large powers of x may then reflect 
an algorithm which pushes the limits of computational 
feasibility. 

The trade-off between sophistication of an Ansatz and 
the associated cost of the update algorithm represents a 
critical tension in the development of novel tensor net- 
work algorithms. For example, the structure of the 4:1 
2D MERA d, nil indicates that this particular Ansatz 
will provide a powerful representation of highly-entangled 



2D systems, and this has been confirmed analytically in 
Ref. |4J where it is shown to furnish a compact and phys- 
ically meaningful description of the toric code. However, 
numerical computations using this Ansatz are hindered 
by an update cost of 0(x^^)- The ability to quickly and 
conveniently determine the cost of contracting different 
tensor networks is therefore of great importance to re- 
searchers employed in the development of novel tensor 
network algorithms. 

Even when the cost of updating a tensor network algo- 
rithm is known, the implementation of these algorithms 
is frequently a non-trivial affair. The process of contract- 
ing a tensor network may always be viewed as a series of 
pairwise contractions, and the overall cost of contracting 
the network is highly dependent upon the sequence in 
which these contractions are carried out. For instance, 
the network shown in Fig.[IJi) is one network which must 
be contracted during the variational optimisation of the 
3:1 ID MERA [Hlgil. The most efficient contraction 
sequence yields a cost of O(x^), but careless choices of 
sequence can yield costs as high as 0(x^^). Thus the 
ability to determine the optimal contraction sequence for 
a tensor network is also important to those implementing 
pre-existing tensor network algorithms, if these are to be 
implemented in a computationally efficient manner. 

Until now, the only option for determining the optimal 
contraction sequence for a tensor network algorithm has 
been a labour-intensive manual study of the tensor net- 
work, and for all but the simplest networks an exhaustive 
search is unfeasible. The process is time-consuming, and 
for networks composed of many tensors with similar num- 
bers of indices but subtly different connectivity, it can be 
difficult to be certain whether any contraction sequence 
is indeed optimal. This paper introduces an algorithm 
which determines an optimal contraction sequence for 
any given tensor network, permitting researchers who de- 
sign and implement tensor network algorithms to quickly 
and easily determine efficient network contractions and 
to identify with confidence the minimal computational 
cost. The algorithm is described in Sec. [Ill and a refer- 
ence implementation is provided in the archive associated 
with this paper and documented in Sec. IIIIl 




II. ALGORITHM 



FIG. 1. (i) Graphical representation of one of the tensor net- 
works which must be contracted during variational optimi- 
sation of the 3:1 MERA [13, ^. (ii) A bipartition of this 
network defines two tensors, X and Y, corresponding to the 
contraction of the resulting sub-networks, as shown. A cost 
may be associated with this bipartition corresponding to the 
cost of contracting X with Y. (iii) The network associated 
with tensor X may then, in turn, be further subdivided into 
networks associated with tensors V and W, as shown, (iv) To 
describe bipartitions of a network, it is useful to label indices. 
In this diagram, summed indices are labelled with unique pos- 
itive integers, and open indices are labelled with unique nega- 
tive integers descending consecutively from —1. Tensors have 
also been labelled with letters to allow easy reference from 
the text. 



non-disjoiniQ and contain no traces fSec. lII At . and then 
allowing disjoint networks and sub-networks fSec. HTbI) . 
Additional techniques used to optimise the search process 
are discussed in Sec. Ill C[ and finally traces are included 
in Sec. HTdI 



This Section describes the Netcon algorithm (short 
for "Network Contraction" ) for performing an optimised 
search through the space of possible tensor network con- 
traction sequences. The problem of determining the most 
efficient contraction sequence for a given tensor network 
is reformulated as a search through the space of sequen- 
tial subdivisions of that network, first under the assump- 
tion that the initial network and all sub-networks are 



^ A tensor network is disjoint iff its tensors can be separated into 
two sub-networks sharing no common indices. A tensor network 
is non-disjoint iff it is not disjoint. 



A. Non-disjoint networks 

1. Basic procedure 

Consider a non-disjoint tensor network such as that 
shown in Fig. [TJi) . This network may be contracted to 
a single tensor by evaluating the sums over paired in- 
dices in some sequence, with the computational cost of 
that sequence corresponding to the number of multiplica- 
tion operations which must be performed. However, not 
all sequences of pairwise contractions attract equivalent 
computational costs. It is shown in Appendix [K\ that for 
any tensor network there always exists a sequence which 
achieves the minimum possible cost using only pairwise 
contractions, and so the objective of the Netcon algo- 
rithm will be the identification of a pairwise contraction 
sequence which realises this minimum cost. 

One might therefore propose an iterative algorithm to 
explore the space of possible contraction sequences, main- 
taining a list of available tensors and their connections 
and updating this list as the tensors undergo pairwise 
contraction. While this approach is viable, the coding 
(and performance) overhead associated with maintain- 
ing both the list of tensors and the history of individual 
tensors, so that composite tensors may be uncontracted 
and different contraction sequences explored, makes this 
approach unattractive. 

The problem of determining the optimal contraction 
sequence for a tensor network such as Fig. [IJi) may be 
recast more favourably as follows: Rather than consid- 
ering pairwise contractions of tensors, instead consider 
bipartitions of the original tensor network. An example 
bipartition is shown in Fig. HJii), and divides the orig- 
inal network (which will be denoted Afz) into two sub- 
networks A/x and TVy (which for now will also be assumed 
to be non-disjoint). The tensors obtained on contracting 
these two sub-networks may be denoted X and Y, and 
the minimum cost of contracting network A/z finishing 
with the contraction of tensor X with tensor Y is then 
the minimum cost of contracting network A/x, plus the 
minimum cost of contracting network A/y, plus the cost 
associated with this final contraction. Provided the min- 
imum costs for contracting all sub-networks such as Ax 
and A/y are known, the minimum cost for contracting 
A/z may then be computed by iterating over all possible 
bipartitions of A/z and identifying the cheapest option. 
Evaluation of the costs to contract sub-networks such as 
Ax and A/y may be performed by recursion as needed, 
with the cost associated with a trivial network (one con- 
taining only one tensor) being zero. If getcost(A/z) is 
a function which performs the above-described iteration 
over bipartitions of A/z, then getcost (A/z) may compute 
the costs associated with contraction of networks A/x and 
A/y by invoking getcost (TVx) and getcost (TVy) respec- 
tively. In this way, a function getcost () may explore all 
possible contraction sequences for a network Az without 
explicitly constructing the dimensions of the intermedi- 
ate tensors. 



In conjunction with this approach, it is useful to intro- 
duce a means of describing bipartitions of a tensor net- 
work which is similarly independent of the intermediate 
tensors involved. For a non-disjoint network partitioned 
into two non-disjoint sub-networks, any such partition 
may be uniquely described in terms of the summed index 
pairs which are split, such that one occurrence of the in- 
dex lies on network A/x and the other lies on network A/y. 
In Fig. [Ijiv) a unique label has been assigned to each in- 
dex of Fig. [T]^i); any convention may be adopted for these 
labels, but in the present paper all summed indices are 
labelled with positive integers and all unsummed indices 
with negative integers descending consecutively from —1. 
The tensors themselves have been labelled with the let- 
ters A. . . G for reference. Given the labelling of Fig.[IJiv), 
contraction of tensor X with tensor Y is associated with 
contraction over the indices labelled 1, 2, 4, 9, 11, and 
12, and this list of indices is independent of how networks 
Ax and A/y might be subdivided in turn. One may then 
subdivide network A/x into two portions, AA/ and A/w, as 
shown in Fig. [ijiii), with contraction of V and W being 
associated with the indices 6 and 7. Given tensors V, W, 
and Y, the index sequence 

6 7 1 2 4 9 11 12 

then specifies that first V is to be contracted with W, 
then the resulting tensor (X) is to be contracted with 
Y. A sequence containing all summed index labels on a 
network A/z [e.g. 1 to 12 in Fig. [Hiv)] thus specifies a 
unique sequence of pairwise contractions whereby that 
network may be fully contracted into a single tensor. 

This method of describing contraction sequences may 
be contrasted with one based on the letter labels of 
Fig. [Tfiv): If brackets () are used to denote pairwise con- 
tractions, so that e.g. (BC) represents the contraction 
of tensor B with tensor C, then tensors X and Y can 
each be represented by any of a number of different se- 
quences depending on the contraction sequence by which 
they are constructed, and the way in which contraction 
of X with Y is represented will vary accordingly. These 
two methods of describing contraction sequences are in- 
terchangeable for non-disjoint networks; for example, for 
Fig- IHiv) the index sequence 

94657213 11 12 8 10 

uniquely specifies the pairwise contraction sequence 

(((((CE)D)A)F)(BG)) 

and vice versa. However, within the Netcon algorithm 
an index-based approach is preferable as the description 
of a bifurcation of network A/z is then independent of 
the optimal contraction sequences of the sub-networks 
A/x and A/y . Once an optimal sequence has been deter- 
mined, this may readily be translated into a tensor-based 
notation if desired (with this process being summarised 
in Appendix [B|) . 
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FIG. 2. Diagram showing part of a tensor network in which 
two indices, labelled 4 and 6, connect the same two tensors. 
Dotted lines represent summed indices which connect the two 
tensors shown here to the rest of the network. 



Note that when two or more consecutive indices in an 
index sequence connect the same tensors, for example 
the indices labelled 4 and 6 in Fig. [51 it is assumed that 
these indices should first be combined into a single index 
having dimension equal to the product of its constituents. 
For example, if the index labelled 4 in Fig. [2] can range in 
value from 1 to 2 (i.e. the index labelled 4 is of dimension 
2) and the index labelled 6 is of dimension 3, these two 
indices may be replaced by a single index t ranging from 
1 to 2 X 3 = 6, as per TableHl Note that this combining of 
indices applies even when composite objects are involved, 
and thus if a contraction sequence for Fig. [Ijiv) begins 
9 4 6 . . . , the first contraction (over index 9) combines 
tensors C and E, generating a new tensor which then 
shares both indices 4 and 6 with tensor D. Indices 4 and 
6 are therefore to be combined into a single index before 
contracting. 

Taking this choice of notation into account, a pseu- 
docode function getcostCAz) which returns both the 
minimum cost associated with contracting a network A/z 
and an index-based contraction sequence realising that 
cost may be represented as follows: 

Define: [sequence cost] =getcost (TVz) 

1: Let the value of cost be marked as unassigned. 

2: Iterating over all bipartitions of TVz, for each bi- 
partition: 
2a: Let A/x and Afy denote the sub-networks 



TABLE I. Example replacement of the indices labelled 4 and 
6 (of dimensions 2 and 3 respectively) with a single index, 
denoted i, of dimension 6. 



Value of index t 


Associated values of indices 4 and 6 


Index 4 


Index 6 


1 
2 
3 
4 
5 
6 


1 
1 
1 
2 
2 
2 


1 
2 
3 
1 
2 
3 



resulting from this bipartition, and X and 
Y denote the tensors obtained on con- 
tracting A/x and A/y • 

2b: Let costl be the cost of contracting ten- 
sor X with tensor Y. 

2c: Let seql be the list of index labels con- 
necting tensor X with tensor Y. 

2c: If A/x contains more than one tensor, 
let [seqX costX]=getcost (A/x), else let 
seqX be empty and let costX=0. 

2d: If A/y contains more than one tensor, 
let [seqY costY]=getcost (A/y), else let 
seqY be empty and let costY=0. 

2e: Let costZ be costl+costX+costY. 

2f: Let seqZ be the sequential concatenation 
[seqX seqY seql] . 

2g: If the value of cost is unassigned, or if a 
value has been assigned to cost and this 
value satisfies cost>costZ, let cost take 
the value of costZ and let sequence take 
the value of seqZ. 
Return the values of sequence and cost. 



2. Iteration over nan-disjoint sub-networks 

Given sufficient time, the algorithm described above 
will always return the minimum cost for contracting any 
tensor network Az and an associated contraction se- 
quence. However, some thought must be given with re- 
spect to the manner in which the algorithm explores the 
space of bipartitions of the network, or else "sufficient 
time" may be very long indeed. In the limit of a fully- 
connected network of n tensors (i.e. one in which every 
tensor is connected to every other tensor), the number 
of bipartitions to be considered in the initial step scales 
as 2"~^ — 1 and the number of admissible contraction 
sequences scales as 



n\{n — 1)! 
2"-i ■ 



(2) 



In practice tensor networks of interest are less intercon- 
nected than this, so restriction to bipartitions yielding 
non-disjoint sub-networks will substantially reduce these 
figures, but the cost of identifying an optimal contraction 
sequence may still be anticipated to rise rapidly with n. 
When the requirement for non-disjoint subnetworks is 
relaxed in Sec. IIIBi the number of possible contraction 
sequences for a tensor network returns to the worst case 
scenario given in Eq. ([2]). However, it turns out that 
only a fraction of these sequences need be taken into ac- 
count in order to be certain of obtaining a sequence of 
minimum cost. It is therefore useful to describe a pro- 
cedure for iterating over all bipartitions of a non-disjoint 
tensor network A/z which yield non-disjoint sub-networks 
Ax and A/y, with this procedure requiring only minimal 
modification in Sec. Ill B 31 



To begin, consider a fully-connected network of n ten- 
sors. For such a network, one may explore all possible 
bipartitions yielding non-disjoint sub-networks using the 
following procedure: First, assign a unique label to each 
tensor (e.g. for n — 7, the letters A. . . G). Write out these 
labels in sequence. Then to iterate over bipartitions, let 
a variable x count from 1 to 2"^^ — 1. For each value of 
X write its binary representation under the list of tensor 
labels. Where a label is associated with a 0, this tensor 
is placed in group one. Where a label is associated with 
a 1, this tensor is placed in group two. This process is 
illustrated in Table HIl for n = 7. The value x = may be 
omitted as for this value group 2 is empty, so this does 
not represent a bipartition of the network. It is also un- 
necessary to consider x G [2'^~^ + 1,2"'] as the resulting 
bipartitions are equivalent to those for x G [0, 2"~^] with 
groups one and two interchanged. 

For a realistic network with fewer interconnections, 
however, following this procedure to bipartition the net- 
work will frequently result in more than two disconnected 
pieces. Iteration over all possible bipartitions as per Ta- 
ble |TT] is therefore relatively inefficient for realistic tensor 
networks, but may be built on to realise a more econom- 
ical iteration over bipartitions of the tensor network as 
follows: 

First, note that if a bipartition such as Group 1 — 
[ABCD], Group 2 = [EFG] represents a valid bipartition 
of the original network (i.e. that neither group is dis- 
joint), then it is unnecessary to consider also Group 1 — 
[EFG], Group 2 = [ABCD]. We may choose to impose 
this constraint by constructing an algorithm which only 
generates bipartitions in which tensor A is in group two. 

Next, define Tier 1 (denoted 7i) as being a set of ten- 
sors whose only member is tensor A, and let tensor A 
be in group two. Now define Tier 2 as the set of all 
tensors in contact with a member of Tier 1 which is in 
group two, but which are not themselves in Tier 1, and 
let all members of Tier 2 initially belong to group one. 
This definition is then extended to higher tiers such that 
Tier x for a; > 1 is the set of all tensors which are con- 
nected to a tensor in Tx-i which belongs to group two, 
but which are not themselves a member of tier Ty for any 
y < X. Initially, Tx will be empty for all x > 2. In gen- 
eral, not all tensors will be assigned to a tier, and tensors 



TABLE II. Using the binary representation of a counter x to 
iterate over all possible bipartitions of a list of seven tensors 
denoted A. . . G. 




X 


A B C D E F G 


Group 1 


Group 2 


1 


1 


ABCDEF 


G 


2 


10 


ABCDEG 


F 


3 


11 


ABODE 


EG 


4 


10 


ABCDFG 


E 




63 


111111 


A 


BCDEFG 



<z> 



FIG. 3. Example tensor network used to demonstrate the 
contiguous-group iteration algorithm described in Sec. Ill A II 
The output of this algorithm is given in Table IIIII 



not assigned to a tier are always taken to belong to group 
one. This initial configuration therefore places tensor A 
in group two, and all other tensors in group one. 

Advancing to the next bipartition configuration is now 
performed as follows: 

1. Let X be the number of the highest non-empty tier 
(initially. Tier 2). 

2. In the manner of Table |lTl perform a binary in- 
crementation of the group assignations in Tier x. 
For example, in Fig. [ijiv) Tier 2 would comprise 
tensors C, D, and F. Initially, all three belong to 
group one. After incrementation, one tensor (say 
F) would now belong to group two. 

3. If incrementing the configuration resulted in binary 
overflow, such that after the increment all elements 
of Tier x belong to group one (e.g. for Tier x con- 
taining three elements, this would occur on incre- 
menting from 111 to 000), then: 

3a. If the tier which overflowed was Tier 1, itera- 
tion is complete. Stop. 

3b. Otherwise, let Tier x be empty, and return to 
Step 1. 

4. Having changed the group assignations in Tier x, 
Tier x -\- 1 may now be non-empty. Let all tensors 
in Tier x -I- 1 be in group one. 

5. Read off new configuration. 

Repeating this process iterates over all bipartitions for 
which group two is non-disjoint. It is still necessary to 
check whether group one is non-empty and non-disjoint, 
but by using this approach, the number of bipartitions 
which must be examined may be substantially reduced. 
A full example is given in Fig. [3] and Table IIIII 



B. Disjoint networks 

Discussion thus far has been restricted to tensor net- 
works which arc non-disjoint. However, sometimes con- 
structing an optimal contraction sequence for a tensor 
network — even one which is non-disjoint — may necessar- 
ily involve taking the product of two tensors which do not 
share any common indices (an "outer product" ) , with an 
example of such a network being seen in Fig. 21 and to 



Iteration 


Tier 1 


Tier 2 


Tier 3 


Tier 4 


No Tier 


Group 1 


Group 2 


Group 1 non-disjoint 
and non-empty? 


1 


Tensors: A 
Groups: 2 


Tensors: B 
Groups: 1 


Tensors: - 
Groups: - 


Tensors: - 
Groups: - 


Tensors: C,D,E 
Groups: 1,1,1 


B,C,D,E 


A 


Yes 


2 


Tensors: A 
Groups: 2 


Tensors: B 
Groups: 2 


Tensors: C,D 
Groups: 1 , 1 


Tensors: - 
Groups: - 


Tensors: E 
Groups: 1 


C,D,E 


A,B 


No 


3 


Tensors: A 
Groups: 2 


Tensors: B 
Groups: 2 


Tensors: C,D 
Groups: 1,2 


Tensors: E 
Groups: 1 


Tensors: - 
Groups: - 


C,E 


A,B,D 


No 


4 


Tensors: A 
Groups: 2 


Tensors: B 
Groups: 2 


Tensors: C,D 
Groups: 1,2 


Tensors: E 
Groups: 2 


Tensors: - 
Groups: - 


C 


A,B,D,E 


Yes 


5 


Tensors: A 
Groups: 2 


Tensors: B 
Groups: 2 


Tensors: C,D 
Groups: 2,1 


Tensors: - 
Groups: - 


Tensors: E 
Groups: 1 


D,E 


A,B,C 


Yes 


6 


Tensors: A 
Groups: 2 


Tensors: B 
Groups: 2 


Tensors: C,D 
Groups: 2,2 


Tensors: E 
Groups: 1 


Tensors: - 
Groups: - 


E 


A,B,C,D 


Yes 


7 


Tensors: A 
Groups: 2 


Tensors: B 
Groups: 2 


Tensors: C,D 
Groups: 2,2 


Tensors: E 
Groups: 2 


Tensors: - 
Groups: - 


- 


A,B,C,D,E 


No 



TABLE III. Iteration over all bipartitions of a tensor network for which group two is not disjoint, using the algorithm described 
in Sec. Ill Cl The tensor network being iterated over is given in Fig. O 



D 



where the original tensor network A/z is itself disjoint. 



Optimal sequences involving outer products 



FIG. 4. An example tensor network for which the optimal 
contraction sequence involves performing an outer product. 
All indices are of dimension x, with x > 1- 



locate such a contraction sequence it is necessary to con- 
sider bipartitions which generate disjoint sub-networks. 
The problems associated with this are two- fold: First, 
it is necessary to extend the index-based contraction se- 
quence notation employed in Sec. Ill A II to encompass 
outer products, and second, it is necessary to identify 
an optimised exploration of the space of contraction se- 
quences which includes disjoint sub-networks when re- 
quired, but which does not involve evaluating the full 
list of 0(n!) possible contraction sequences described in 
Sec HTO 

Initially considering only networks A/z which are non- 
disjoint. Sec. Ill BT] examines the conditions under which 
it is necessary to consider outer products in order to 
be certain of finding a sequence of optimal cost, and 
Sec. Ill B"2] introduces an extension of the index-labelling 
notation of Sec. Ill AT] which is capable of describing the 
resulting contraction sequences. It is shown in Sec. IIIB II 
that even when outer products (and hence bipartitions 
which yield disjoint sub-networks) are taken into account, 
it is still not necessary to examine all bipartitions of a net- 
work TVz, and Sec. IIIB 31 then explains how the iteration 
algorithm of Sec. Ill A 21 may be modified to include the 
necessary additional bipartitions without needing to it- 
erate over all bipartitions of a non-disjoint network A/z. 
Finally, Sec. Ill B 41 extends this approach to situations 



As mentioned in Sec. Ill A 21 the performance of any 
algorithm which attempts to iterate over all bipartitions 
of a tensor network, including bipartitions yielding dis- 
joint sub-networks, will scale very poorly with the num- 
ber of tensors involved. Fortunately, it is not necessary 
to consider all bipartitions — only those which are neces- 
sary in order to be certain of finding an optimal contrac- 
tion sequence — and it is consequently possible to derive 
a number of restrictions on which bipartitions yielding 
disjoint subnetworks should be taken into consideration. 

First, note that performing an outer product between 
two tensors is formally equivalent to contracting over a 
shared index of dimension 1 . Now consider a contraction 
sequence in which an outer product is taken between two 
tensors, A and B, and then the resulting object (AB) is 
contracted with a further tensor C which shares one or 
more indices with A but none with B. This situation is 
illustrated in Fig. [S] in which 

• all indices on tensor A not connecting to tensor C 
have been combined into a single index of dimen- 
sion a, 

• all indices on tensor C not connecting to tensor A 
have been combined into a single index of dimen- 
sion c, 

• all indices connecting tensors A and C have been 
combined into a single index of dimension d, 

• and all indices on tensor B have been combined into 
a single index of dimension b. 

As the combination procedure (which is illustrated in Ta- 
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FIG. 5. An example sub-network whose contraction involves 
an outer product. Should tensor A be contracted with ten- 
sor C before or after performing an outer product involving 
tensor B? Letters a, . . . ,d denote index dimensions. 
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FIG. 6. Two tensors A and B share no common index; should 
they be combined using an outer product before contracting 
with tensor C? Letters a, . . . , e denote index dimensions. 



ble|T| is reversible, this may be done without loss of gen- 
erality. Tensors A and B share no common indices by 
definition, as the contraction of A with B is required 
to be an outer product. Given this labelling, one can 
ask under what circumstances the contraction sequence 
((AB)C) may be superior to the sequence ((AC)B). Eval- 
uating the costs of these two contraction sequences yields 
the inequality 



abd + abed < acd + abc 



(3) 



which, for a, 6, c, and d being positive integers, is satisfi- 
able only if fe = 1 or rf = 1 . Defining the total dimension 
of a tensor as the product of the dimensions of all its 
indices, it follows that performing an outer product be- 
tween any two tensors A and B may only ever potentially 
be advantageous if: 

• the next contraction to be performed is also an 
outer product (equivalent to d = 1), 

• the next contraction is with a tensor C that shares 
non-trivial indices (i.e. indices of dimension greater 
than 1) with both A and B, or 

• one or more of tensors A and B has total dimen- 
sion 1 (i.e. is a single number). 

Networks involving tensors of total dimension 1 repre- 
sent a special case of limited interest, and may always be 
factored out as a simple numerical multiplier. It is conve- 
nient to neglect these cases in the interest of obtaining a 
better-performing algorithm for tensor networks of prac- 
tical interest. 

Next, consider the most general situation in which ten- 
sor C shares non-trivial indices with both A and B. This 
situation may be represented as shown in Fig. [51 where 
d > 2 and e > 2. All indices on tensor A not connect- 
ing to tensor C have been combined into a single index 
of dimension a, and similarly for b and c. All indices 
common to tensors A and C have been combined into a 
single index d, and all those common to B and C have 
been combined into a single index e. As before, tensors 



A and B share no common indices by definition. The 
three possible contraction sequences for this diagram and 
their associated costs are now given in Table HVl and re- 
quiring that ((AB)C) be cheaper than both of the other 
two (permitting an outer product to potentially offer an 
advantageous contraction sequence) gives rise to the in- 
equalities 



c{b + d- bd) > bd 

c{a + e — ae) > ae. 



(4) 

(5) 



Since all indices must take positive integer values greater 
than or equal to one, and d and e are required to be 
greater than or equal to two, it is seen that these equa- 
tions admit solutions only for 



1, 



reducing to 



c> d 
c > e. 



(6) 



(7) 
(8) 



Consequently, if the optimal contraction of tensors A, B, 
and C must begin with an outer product between A and 
B, tensors A and B may not have any non-trivial indices 
other than those which they share with tensor C, and 
tensor C must have at least one such index, of dimension 
satisfying Eqs. ^ and (0). 

In combination with the results inferred from Fig. [S] 
it follows that sequences involving an outer product on 
n tensors Ai , . . . , A„ need only be considered when they 
satisfy the following criteria: 

• The outer products are performed pairwise (see Ap- 
pendix El . 

• After performing n — 2 pairwise outer products, the 
series of outer products yields a pair of tensors A 
and B each having one or more non-trivial indices 
shared with tensor C and no non-trivial indices not 
shared with tensor C. 

As it is never optimal to perform an outer product be- 
tween two tensors which share a non-trivial summed in- 
dex (as opposed to summing over that index), and as 
networks involving tensors of trivial total dimension have 
been excluded, it follows that all the tensors Ai, . . . , A„ 
each have one or more non-trivial indices shared with 
tensor C and no non-trivial indices not shared with ten- 
sor C. 



TABLE IV. Contraction sequences and associated costs for 
the tensor network shown in Fig. [6] 



Sequence 


Cost 


((AB)C) 
((AC)B) 
((BC)A) 


abde + abcde 
acde + abce 
bcde + abed 
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Next, consider the construction of tensor C. It turns 
out that Eqs. dl])-® permit us to conclude that a con- 
traction sequence for tensor C can always be chosen such 
that it is not itself the result of an outer product, while 
still obtaining an optimal cost overall; If 

• the optimal final contraction in the construction of 
C is unavoidably an outer product of two tensors 
A and B, 

• the optimal contraction of A, B, and C begins with 
the outer product of A and B, and 

• the product of contracting A with B is denoted C, 

then in addition to Eqs. ([5])-® the indices of tensors A, 
B, and C must also satisfy the constraints of Eqs. ©-([S]), 
implying 



c=l 

ab > d 
ab > e. 



(9) 
(10) 

(11) 



Simultaneous satisfaction of Eqs. (H])-® and Eqs. ([9])- 
(llip is impossible, implying that in the construction of an 
optimal contraction sequence it is never necessary to con- 
tract together two tensors both produced using an outer 
product. Consequently, when exploring bipartitions of a 
tensor network A/z it is still permissible to impose that 
one of the two groups of tensors generated by a bipar- 
tition be non-disjoint. If the other is disjoint, it is fur- 
ther permissible to discard this bipartition if the disjoint 
group possesses any unpaired indices other than those 
which are to be contracted with the non-disjoint group. 
Finally, when a sequence of outer products is to be 
performed on a set of tensors Ai,...,A„, the optimal 
order of contraction may be determined by associating 
with each tensor a total dimension ^i, and proceeding as 
follows: 

1. Identify the two tensors having the smallest total 
dimensions. 

2. Remove those two tensors from the list and replace 
with their outer product. 

3. Repeat until only one tensor remains. 

For example, having four tensors Ai , . . . , A4 of total di- 
mensions ^1 = x'': 6 = X^: 6 = X*^, ^4 = X^ for some 
value Xy 3-^1 optimal sequence for taking the outer product 
of these tensors is ((AiA2)(A3A4)). The resulting tensor 
is then contracted with tensor C. 



2. Notation 

Addressing first the problem of notation, consider 
again the example given in Fig. Ul Using a tensor-based 
notation, for this example the optimal contraction se- 
quence may be represented 



Tensors (AB) and C share no common indices, and thus 
the contraction of (AB) with C is an outer product. In 
index-based notation the order over which the labelled 
indices are to be contracted may be written 



12 3 4 5 



(13) 



((((AB)C)D)E). 



(12) 



This sequence, however, makes no reference to the outer 
product which must be performed, and attempting to 
reconstruct Eq. P^ will erroneously yield the subop- 
timal contraction sequence ((((AB)D)C)E). Admittedly 
the sequence of Eq. ([T^ is not unique; recall that 
when multiple indices connect the same two tensors, 
they will be combined into a single index before con- 
traction and hence their ordering is immaterial. How- 
ever, the alternatives fare no better. Like Eq. P^ 
the index sequence 12 3 5 4 implies a tensor con- 
traction sequence ((((AB)D)C)E), and the other two 
alternatives — 1 3 2 4 5 and 13 2 5 4 — fare no bet- 
ter, implying the sequence (((AB)(CD))E) which is again 
sub-optimal. Some extension of index-based sequence no- 
tation is clearly required. 

In Sec. Ill BT] it was shown that to be certain of identi- 
fying a contraction sequence of minimum cost, one need 
only take into account sequences where all tensors A, B, 
etc. participating in an outer product share a summed 
index with some other tensor X which does not partic- 
ipate in the outer product, and the next operation is 
to contract the result of the outer product with X. The 
presence of an outer product of n tensors may therefore 
be denoted by inserting n — 1 zeros into the contraction 
sequence. Conversely, when interpreting a contraction 
sequence, on encountering a string of n — 1 zeros the re- 
sponse is then as follows: 

1. Read the next x indices of the contraction sequence, 
where x is the maximum value such that n tensors 
will be involved in summing over all x indices. 

2. Identify the n tensors involved. 

3. One of these tensors will be connected to all the 
others; this is tensor X as described above. 

4. Perform an optimal pairwise outer product on all 
tensors except X. 

5. Contract the resulting object with tensor X. 

This process is usefully supplemented by two additional 
sub-algorithms. The first is the process described in 
Sec. IIIB II for determining the optimal ordering of a se- 
quence of outer products, and the second identifies X: 

1. Consider the first of the x indices. This index will 
connect two tensors, to be denoted A and B. 

2. Advance through the x indices until one is found 
which appears on A or B, but not both. The ten- 
sor (either A or B) on which this index appears is 
tensor X. 

Provided index-based sequence notation is only used 
to describe outer products of the sort discussed in 
Sec. IIIB II it follows that any contraction sequence con- 
sistent with these restrictions may be described using the 



zero notation described above. The bipartitioning proce- 
dure employed by the Netcon algorithm may be chosen 
so as to yield all contraction sequences which are com- 
pliant with these restrictions, and only contraction se- 
quences which are compliant with these restrictions (see 
Sec lIIBB]) . and since there always exists an optimal con- 
traction sequence which is consistent with these restric- 
tions (see Sec. IIIBTj) . it follows that Netcon will always 
find an optimal sequence, and that the optimal contrac- 
tion sequence returned by Netcon may always be de- 
scribed using the zero notation. This notation is slightly 
less intuitive than the pairwise bracketing notation used 
in (for example) Eq. ([12]), but has the advantage of be- 
ing a linear rather than heirarchically-structured repre- 
sentation of the contraction sequence. For easy reference 
the interpretation of an index-based contraction sequence 
which may include the zero notation is summarised in 
Appendix [B] 

For the example given in Fig. 21 the index se- 
quence encoding the optimal tensor contraction sequence 
((((AB)C)D)E) is 

10 2 3 4 5. 



Completion of the iteration algorithm of Sec. Ill A 21 
amended as above, corresponds to iterating over: 

• all bipartitions for which both groups are non- 
disjoint, and 

• all bipartitions in which the group containing ten- 
sor A is non-disjoint but the other group is disjoint. 

Next, move tensor A from Tier 1 to Tier and select a 
different tensor (e.g. B) to form a new Tier 1. Repeat the 
iteration algorithm of Sec. Ill A2[ but this time examine 
only bipartitions for which group one is disjoint (parti- 
tions for which both groups are non-disjoint have already 
been examined during the iteration in which Tier was 
empty) . On completion the contents of Tier 1 are again 
added to Tier and a new Tier 1 is selected. This pro- 
cedure repeats until all tensors are in Tier 0. These ad- 
ditional steps yield 

• all bipartitions in which the group containing ten- 
sor A is disjoint but the other group is non-disjoint, 

completing iteration through all necessary bipartitions of 

ATz. 



3. Iteration over disjoint and non-disjoint sub-networks of 
a non-disjoint network 



Given the results of Sec. IIIB II when searching for the 
optimal contraction sequence for a non-disjoint network 
A/z it is only necessary to iterate over all bipartitions such 
that at least one of the resulting groups is non-disjoint. 
This may be achieved as follows: 

Let there be an additional Tier not discussed in 
Sec. Ill A^ named Tier 0, and initially empty. All tensors 
in Tier are always in group one. Now select a tensor A 
to use as Tier 1, and perform the iteration algorithm de- 
scribed in Sec. Ill A 21 but do not discard bipartitions in 
which group one is disjoint unless it also possesses open 
indices not shared with group two. Let A/x denote the 
sub-network which is allowed to be disjoint, and replace 
step 2c of the iteration algorithm with 

2c: If A/x contains only one tensor, let seqX be empty 
and let costX=0, otherwise: 

2cl: Let n be the number of disjoint sections mak- 
ing up sub-network A/x, and let A/xi ,. ■ • i A/x„ 
enumerate these sections. 

2c2: Let costOP be the cost for taking the outer 
product of these disjoint sections according to 
the method of Sec. IIIB II 

2c3: For each disjoint section A/xi, let 
[seqXi costXj] =getcost (Axi) . 

2c4: Let seqX be the concatenation 
[seqXi seqX2 . . . seqX„] . 

2c5: Let costX be costOP -I- J^i costX^. 

2c6: Prepend n—1 zeros onto seql. 



4- Computing the optimal contraction sequence for a 
disjoint network 

Up to this point, it has been assumed that the origi- 
nal tensor network A/z supplied to the Netcon algorithm 
is non-disjoint. While iterating through bipartitions of 
A/z the function getcostO may yield disconnected sub- 
networks Ax, but each of these is then further separated 
into its non-disjoint components and getcostO is recur- 
sively applied to each non-disjoint component in turn. 
Extension to situations where the original network A/z is 
disjoint is not difficult. 

Consider again the argument presented in Sec. IIIB II 
regarding the timing of outer products. As it is always 
possible to defer an outer product until such a time as 
the next contraction involves indices shared with both el- 
ements of the outer product, it follows that for a disjoint 
network A/z, taking the outer product between the vari- 
ous elements of A/z should be the final step in contracting 
the tensor network. The process of determining an opti- 
mal contraction for a network A/z composed of n disjoint 
sub-networks may therefore be broken down into that 
of determining an optimal contraction sequence for each 
disjoint sub-network, concatenating these sequences, and 
then appending n — 1 zeros onto the end of the sequence 
to indicate the performance of n — 1 outer products be- 
tween the remaining n tensors. The interpretation of a 
string of zeros in the contraction sequence is therefore as 
follows: 

• When a string of n — 1 zeros is encountered, this 
indicates that n—1 outer products are to be per- 
formed. 



10 



• If only n tensors remain, then the outer product 
of all n tensors is to be evaluated pairwise in the 
optimal order according to the procedure given in 
Sec HTBTl 

• If more than n tensors remain, the n tensors to be 
contracted together are to be determined according 
to the algorithm given in Sec. Ill B 21 as is tensor X. 
The outer product of these n tensors is then to be 
evaluated pairwise in the optimal order according 
to the procedure given in Sec. IIIB II 

Note that some caution is required if the tensor net- 
work Afz includes indices of dimension 1, since (as noted 
in Sec. IIIB 1[) contracting over a shared index of dimen- 
sion 1 is formally equivalent to performing an outer prod- 
uct between two tensors, and should be treated as such. 
This is particularly important where a tensor network 
A/z is disjoint, for example being made up of two disjoint 
sub- networks Nx and JVy, and one of these sub-networks, 
say A/x, is made up of a further pair of sub- networks, Afv 
and A/w , connected only by indices of dimension 1 . Writ- 
ing V, W, X, and Y for the tensors formed on contracting 
sub- networks A/y, A/w, A/x, and A/y respectively, if sub- 
networks A/x and A/y are each contracted independently 
and then an outer product is performed between tensors 
X and Y, an optimal sequence for the contraction of sub- 
network Ax may be identified in which the final pairwise 
product is an outer product between tensors V and W. 
The final steps in the contraction of network A/z are then 
necessarily ((VW)Y). However, it is possible that the op- 
timal sequence for contracting tensors V, W, and Y is in 
fact either ((VY)W) or ((WY)V), and thus the optimal 
sequence for network Az may be missed. To ensure iden- 
tification of the optimal contraction sequence for Az it 
is necessary to identify all sub-networks which are effec- 
tively disjoint, in this case TVy, A/w, and A/y, indepen- 
dently contract each of these sub-networks, then deter- 
mine the optimal contraction sequence for the resulting 
objects (which are permitted to be of trivial total di- 
mension) according to the algorithm given in Sec. IIIB II 
The identification of all effectively disjoint sub-networks 
is easily achieved by stripping from the network all in- 
dices of dimension 1. An optimal contraction sequence 
for the original network may then be achieved by com- 
puting an optimal sequence for the reduced network and 
appending to this sequence all labels which correspond 
to indices of trivial dimension. 

In this approach, contraction over trivial indices is de- 
ferred until contraction of network A/z is otherwise com- 
plete, at which point they represent traces over indices 
of dimension 1. Performing such a trace incurs zero 
computational cost. In general it is unwise to leave in- 
dices uncontracted — if a pair of tensors A and B share 
a common index of dimension d and that index is left 
unsummed when evaluating (AB), the effect is to unnec- 
essarily increase the total dimension of (AB) by a factor 
of (P, and similarly increase the cost of any subsequent 
contractions involving tensor (AB). When d = 1, how- 
ever, no additional cost is associated with this deferral. 




Index label 


Dimension 


1 


6X 


2 


X 


3 


2 

X 


-1 


6X 


-2 


X 


-3 


2 

X 



FIG. 7. A simple network for which the optimal contraction 
sequence varies with the value of x- 



It is further noted that when A/z is non-disjoint, strip- 
ping of trivial indices in this manner may cause the net- 
work to separate into two or more disjoint components. 
Although this process is not obligatory for non-disjoint 
A/z if the Netcon algorithm is to return an optimal con- 
traction sequence it is harmless, and if the stripping 
causes Az to separate into disjoint components then the 
ability to exploit this disjoint structure yields a substan- 
tial performance increase for the Netcon algorithm. Con- 
sequently, the default behaviour of the Netcon algorithm 
is to ignore the existence (and thus defer contraction) of 
trivial indices even for non-disjoint networks A/z . 



C. Efficient exploration of network bipartitions 

1. Excluding high-cost contraction sequences 

When exploring all possible bipartitions of a network 
Az , a substantial performance increase may be obtained 
by excluding those leading to higher-cost contraction se- 
quences. How this is best done depends upon the manner 
in which the dimensions of the tensor indices are to be 
specified, with two different approaches being considered 
in this paper: 

1. In the first approach, index dimensions are spec- 
ified as integer values. The returned contraction 
sequence then minimises the cost for these specific 
values. 

2. In the second approach, index dimensions are speci- 
fied in the form ax'' for some parameter x, and the 
sequence returned is one which is optimal in the 
limit of large x- A cost of the form ax'''^^ is there- 
fore always larger than a cost of a'x^ regardless of 
the values taken by a and a'. 

The contraction sequences returned by these approaches 
will not necessarily always coincide. For example, con- 
sider the simple network shown in Fig. [71 Contracting 
this network according to the index sequence 3 2 1 in- 
curs a cost of 72x^, while the sequence 12 3 incurs a 
cost of 12x^. The sequence 3 2 1 is therefore clearly to 
be preferred in the limit of large x, but sequence 1 2 3 is 
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preferred for explicit index costs corresponding to values 
of X < 6. 

Consider first the task of excluding higher-cost con- 
traction sequences when index dimensions are specified 
as integers. Let ^cap denote the upper bound on con- 
traction sequences to be explored, and re-examine the 
getcostCTVz) algorithm described in Sec. Ill A II Instep 
2, the algorithm iterates over bipartitions of network TVz 
into sub- networks TVx and A/y, whose contractions yield 
tensors X and Y. The minimal contraction cost associated 
with each such bipartition is then the cost of combining 
tensors X and Y, costl, added to the minimal costs as- 
sociated with contracting each of the sub-networks A/x 
and A/y (costs costX and costY respectively). Evalua- 
tion of these latter costs proceeds by recursion and may 
be lengthy if the number of tensors in either of these sub- 
networks is large. However, if the value of costl exceeds 
^cap then it is unnecessary to proceed with the evaluation 
of either costX or costY as the minimal cost associated 
with this bipartition, costl + costX + costY, also neces- 
sarily exceeds ^cap- Similarly, if costl is acceptable but 
costl + costX exceeds ^cap then evaluation of costY is 
unnecessary. 

When evaluating the minimal cost to contract a tensor 
network Az, it is entirely possible that no contraction se- 
quence will be found which is compatible with the upper 
bound on cost imposed by fcap- Furthermore, this situa- 
tion may be encountered at any level of recursion (for ex- 
ample, the cost associated with a particular bipartition of 
TVz may be acceptable, but there may not exist any con- 
traction sequences for TVx which are cheaper than ^cap)- 
In this situation the getcostCTVz) algorithm should re- 
turn a flag indicating failure to find an acceptable con- 
traction sequence, and also the smallest cost for which 
further investigation of a contraction sequence was dis- 
carded. This value will be denoted ^faii and corresponds 
to the smallest of the following: 

1. The smallest value of costl which exceeded fcap- 

2. The smallest value of costl+f ailX where f ailX is 
the returned value of ^faii for sub- network TVx- 

3. The smallest value of costl+costX which exceeded 

Scap- 

4. The smallest value of costl+costX+f ailY where 
failY is the returned value of ^faii for sub-network 
AV. 

5. The smallest value of costl+costX+costY which 
exceeded ^cap- 

This value, ffaii, then corresponds to the minimum value 
to which ^cap niust be increased in order to permit fur- 
ther evaluation of some of the discarded contraction se- 
quences. Conversely it is also possible that a contraction 
sequence will be found for a cost less than ^cap, in which 
situation ^cap should be decreased accordingly as more 
expensive sequences need no longer be considered. 

Next, consider the task of excluding higher-cost con- 
traction sequences when index dimensions are specified 
in the form ax^- For large Xj ^ cost scaling as 0(x") 



will always be greater than one scaling as 0(x"~^), and 
a simple means of placing an upper bound on the cost 
of contraction sequences to be investigated is to impose 
an upper limit on admissible powers of x- Since the to- 
tal cost of a sequence is just the sum of the costs of its 
constitutent pairwise contractions, this constraint may 
in turn be imposed by placing an upper limit ^powcrcap 
on the maximum acceptable power of x to appear in the 
cost of a single pairwise tensor contraction. (Also note 
that individual pairwise contractions are always mono- 
mial in X-) If the cost costl associated with a given 
bipartition of a tensor network Az into sub-networks A/x 
and TVy involves a power of x greater than ^pcap then the 
costs associated with sub-networks Ax and Ay need not 
be evaluated, and getcostO can advance to the next 
candidate bipartition. 

A similar condition again applies to the cost associated 
with contracting sub-network A/x, such that if costX 
contains a power of x greater than some value ^pcap 
then evaluation of the cost to contract sub-network A/y 
is unnecessary. However, in contrast with the situation 
where index dimensions are specified as integer values, 
it is no longer necessary to explicitly check this also for 
costl -|- costX as if costl is of order at most 0(x^'"''''') 
and costX is of order at most ©(x^'"'"'") then the same is 
necessarily true of costl -|- costX. In the event that no 
contraction sequence exists which has a cost of at most 
Q^'-^Cpcap')^ getcostO should return a value ^faii indicat- 
ing the minimum value to which ^pcap must be increased 
to permit the further evaluation of some discarded con- 
traction sequences. In this approach, the value of ffaii is 
given by the smallest of the following: 



1. The power of x appearing in the smallest value of 
costl which exceeds ©(x^^""''). 

2. The smallest value of ^faii returned by a call to 
getcostCAx)- 

3. The smallest value of ^faii returned by a call to 
getcost (A/y) . 



Finally, if a sequence is identified having a cost of O(x^) 
where ^ < ^pcap then the value of ^pcap should be de- 
creased accordingly. 

Regardless of whether index dimensions are specified as 
integers or in the form ax^, imposing an appropriate cap 
on the contraction sequences to be evaluated can result 
in a significant performance increase. However, choosing 
an appropriate value for ^cap or ^pcap rnay be problem- 
atic if the cost of the optimal contraction sequence is not 
known in advance. Setting this value too high may re- 
sult in suboptimal performance of the Netcon algorithm, 
and setting it too low will result in failure to identify 
an acceptable contraction sequence. A resolution to this 
problem is presented in Sec. IIIC 21 
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2. Caching of sub-network outcomes 

During the recursive bipartitioning of a tensor net- 
work A/z described in Sec. Ill All the function getcostO 
wiU frequently be called upon to evaluate a sub-network 
which has been encountered before. [For example, one 
might first partition Fig. [IJiv) into A/abcfg and A/de, 
or one might split it into A/bcg and Aadef then split 
the latter into JVaf and A/de- Both approaches in- 
volve evaluating the cost associated with contracting sub- 
network A/de-] a significant performance increase may 
be obtained by caching results for previously-solved sub- 
networks, such that on the second and later encounters 
their optimal cost and a sequence yielding this cost may 
be looked up rather than being calculated by invoking 
getcostO for the sub-network. 

Even greater advantage may be obtained when the 
caching of costs for sub-networks is combined with the ex- 
clusion of high-cost contraction sequences as described in 
Sec. IIIC 11 If the value of ^cap (or Cpcap, as appropriate) 
is initially set to a value known to be too low to permit 
successful contraction of the network, getcostO will fail 
and will return the lowest-valued cost for which further 
evaluation was aborted. If ^cap is then increased to this 
value, getcostO is called again, and the cached results 
from the last attempt are retained, with this procedure 
being repeated until a successful contraction sequence is 
found, then the net effect is to automatically exclude 
from evaluation all sequences involving a single contrac- 
tion cost greater than some threshold ^cutoff, where ^cutoff 
corresponds to the largest single contraction cost arising 
during the optimal contraction sequence. Where more 
than one optimal contraction sequence exists, let i enu- 
merate these sequences and let ^^^x be the maximum 
cost associated with a single contraction in sequence i. 
The value of ^cutoff is then the smallest value in the set 

ir }. 

LSinaxJ 

In addition to caching the costs and sequences for 
successfully-evaluated sub-networks, when employing a 
cutoff ^cap it is also useful to cache the costs at which eval- 
uation for a given sub-network is aborted (^faii)- Thus if 
evaluation of sub-network TVp is aborted due to a min- 
imum step cost of at least p > ^cap (i-e. the cheapest 
decomposition of A/p has been determined to inevitably 
involve a contraction having a cost of at least p) and sub- 
network TVq is similarly aborted due to a minimum cost 
q > p, then if fcap is increased to p it is immediately 
apparent that sub-network Ap should be re-evaluated as 
it may perhaps now be contractible with all steps cost- 
ing less than ^cap, whereas there is no need yet to fur- 
ther investigate sub- network A/q. If re-evaluation of sub- 
network A/p shows that it can indeed be contracted with 
a maximum cost-per-step of ^cap then the corresponding 
cost and sequence are added to the cache. If not — if it 
turns out that after performing the step of cost p a more 
expensive contraction is also required — then the cached 
value of ^faii for Ap should be updated instead. 

Using this approach, it is interesting to compare the or- 



der with which previously-excluded sub-networks become 
available for further investigation when index dimensions 
are specified as integers, versus when they are specified 
in the form ax''. For instance, consider a situation where 
index dimensions are specified in the form ax'^ and initial 
investigation of two sub-networks A/p and Aq shows that 
contraction of Ap will involve at least one step costing 
X^ , while contraction of A/q will involve at least one step 
costing 2x^- Under this approach, both sub- networks 
will become eligible for further investigation at the same 
time (when ^pcap is increased to 6) . Whether Ap or A/q 
is evaluated first will depend upon which is encountered 
first when iterating through bipartitions of the parent 
network. If, however, an explicit value is substituted for 
X, say X = 5, and the Netcon algorithm is invoked for 
the same network but with the index dimensions speci- 
fied as integers, the increase of ^cap is such that bipar- 
titions and sub-networks arc investigated in strict order 
of cost. Sub-network A/p will be investigated first, when 
^cap rises to 5^, and A/q will only be investigated if no 
solution exists for a cost less than 2x5^, with this taking 
place on a later iteration of the Netcon algorithm once 
Ccap has been raised to 2 x 5^. In practice, the former 
approach strikes a better balance between making candi- 
date sequences available more rapidly (and so reducing 
the number of iterations of the Netcon algorithm before 
a solution is found) and the exclusion of large numbers of 
unnecessarily costly sequences which may cause individ- 
ual iterations of the Netcon algorithm to take an exces- 
sively long time. Consequently, where it is appropriate 
to do so, specification of index dimensions in the form 
ax'' is to be preferred. In the reference implementation 
provided with this paper, the slower performance when 
index dimensions are specified as integers has been offset 
by requiring that any time ^cap is incremented, the new 
value is at least equal to the old value of ^cap multiplied 
by the dimension of the smallest non-trivial index of the 
tensor network. 



D. Traces 

It is worth making brief mention of tensor networks in- 
volving traces. Unlike tensor contractions, which require 
multiplication, traces may be evaluated using only addi- 
tion. As this operation is relatively cheap, it is assumed 
that the costs associated with tracing are negligible when 
compared to those associated with pairwise contraction. 
Consequently, the Netcon algorithm performs all traces 
on a network before any pairwise contractions. Further, 
as noted in Sec. Ill B 41 when two tensors share multiple 
indices of non-trivial dimension it is never advantageous 
to leave one or more of these indices uncontracted while 
combining the tensors. Consequently, the only traces 
ever appearing in the optimal sequences returned by the 
Netcon algorithm are those where the initial network con- 
tains traces, which are performed at the beginning of the 
contraction sequence, and those over indices of dimen- 
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sion 1, which are performed at the end of the contraction 
sequence as discussed in Sec. Ill B 41 



III. REFERENCE IMPLEMENTATION 

An example implementation of the Netcon algorithm 
described in Sec. HIl written in Matlab and C++, 
may be obtained as follows: While viewing the abstract 
page for this paper on arXiv.org, select "Download > 
Other formats", then "Download source". On expand- 
ing the resulting archive, the files comprising an imple- 
mentation of the Netcon algorithm are netcon. m and 
netcon_cpp. cpp. The implementation requires Matlab 
2011a or above and a compatible C++ compiler, and has 
been tested imder Matlab 2011a with both Microsoft 
Visual C++ 2010 and Apple XCode 4.5.2. 



A. Compilation 

The reference implementation of the Netcon algorithm 
comprises a Matlab function netconO which may be 
invoked from the Matlab command line. It requires 
compilation of a C-l — h component before use, which is 
achieved using the Matlab mex command. First, en- 
sure that a supported C++ compiler has been config- 
ured using the command mex -setup as described in the 
Matlab documentation. Then download and extract 
the contents of netcon.zip to obtain the files netcon. m 
and netcon_cpp. cpp. The C-\ — h component of the im- 
plementation may then be compiled using the command 

netcon compile 

which is equivalent to invoking mex netcon_cpp.cpp. 

Note that by default the C++ portion of the code 
makes use of automatic memory allocation and dealloca- 
tion for variable-length arrays, and that this feature is not 
supported by some C-l — h compilers, including Microsoft 
Visual C++ 2010. If compilation using the default con- 
figuration fails, an alternative configuration which em- 
ploys manual memory allocation and deallocation may 
be compiled using the command 

netcon compile safe 

(equivalent to mex -DSAFE netcon_cpp.cpp). This con- 
figuration, however, is associated on some systems with 
a performance decrease by a factor of 2-3, and thus when 
supported the former configuration is preferred. 



B. Invocation 

Invocation of the algorithm is via the Matlab com- 
mand 

[sequence cost] = netconClegLinks .verbosity, 
costType,xiCap,legCosts .keepTriv) ; 



and takes between one and six input parameters, as fol- 
lows: 

legLinks: This parameter describes the tensor net- 
work for which an optimal contraction sequence is sought. 
To construct legLinks, first draw the tensor network 
using the customary graphical notation (summarised in 
§1.2 of Ref . I45I) . with each tensor being represented by a 
shape, each summed index by a line connecting the two 
tensors on which it appears, and each unsummed index 
by a line with one free end and the other end attached to 
the tensor on which it appears [e.g. Fig.[IJi)]. Next, label 
each summed index with a unique positive integer, and 
each unsummed index with a unique negative integer, de- 
scending consecutively from —1 [e.g. Fig.[IJiv)]. For each 
tensor T, now construct a 1 x ut matrix where nx is the 
number of indices attached to tensor T, with entries cor- 
responding to the labels associated with those indices. 
Ordering of the indices is unimportant. If there are m 
tensors, then there are m such matrices, which may be 
denoted Mi, i G {1 . . .m}. Finally, legLinks comprises 
a 1 X m cell array, with the m entries of this array cor- 
responding to the matrices Mi, with ordering once again 
being unimportant. For example, the labelling given in 
Fig. [ijiv) may be associated with the input parameter 

legLinks = {[-1 1 2 3] , [2 4 5 6] , [1 5 7 -3] , 

[3 8 4 9] , [6 9 7 10] , [-2 8 11 12] , 
[10 11 12 -4]}. 

verbosity: Determines the level of output generated 
by netconO. For verbosity — 0, an optimal index 
contraction sequence is returned in sequence and the 
associated cost is returned in cost but operation is oth- 
erwise silent. For verbosity = 1, a message is gener- 
ated every time the upper bound on tensor contraction 
costs is increased as described in Sec. IIIC 1[ and on com- 
pletion an optimal sequence and the associated cost are 
displayed on screen. For verbosity = 2, behaviour is as 
for verbosity — 1 but candidate contraction sequences 
and associated costs are announced when these sequences 
constitute the lowest-cost contraction sequence found so 
far. The final sequence and cost announced then corre- 
spond to the optimal solution reported at verbosity = 1 
and returned in the output variables sequence and cost. 
Default value: 2. 

costType: To determine an optimal contraction se- 
quence associated with a given tensor network, it is nec- 
essary to specify the dimension of each index in the net- 
work. Index dimensions may either be specified as inte- 
gers, or in the form ax , where x is an unspecified param- 
eter which is presumed to be large. The former is indi- 
cated by costType = 1, and the latter by costType — 2. 
Default value: 2. Note that when using costType — 2, 
indices of fixed dimension d (such as the physical indices 
of an MPS or PEPS) may be represented by setting a = d 
and 6 = 0. 

xiCap: When searching for an optimal contraction se- 
quence, netconO initially restricts itself to sequences 
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having a cost of at most ^cap (for costType = 1) or 
0(x^p''''p) (for costType = 2), where xiCap represents 
Ccap or ^pcap respectively. The value of xiCap will auto- 
matically increase if no contraction sequence exists which 
satisfies this constraint. Setting xiCap too high can incur 
extremely large overheads, whereas the process of auto- 
matic increase is relatively low-cost due to the caching of 
data described in Sec. lIIC2l It is therefore recommended 
that xiCap be left at its default value of 1 unless the cost 
to contract the network is already known. 

legCosts: The format of this input parameter is de- 
pendent upon the value of costType. 

1. For costType = 1, legCosts is an ^ x 2 matrix 
whose first column consists of index labels and 
whose second column gives the dimensions associ- 
ated with those labels. If legCosts is not specified, 
it is assumed that each index has dimension 2. 

2. For costType — 2, legCosts comprises an £ x 3 
matrix where £ is the total number of unique in- 
dex labels. Each row then comprises three entries, 
[x o 6], where a; is a index label (and hence a pos- 
itive or negative integer), and a and b specify the 
dimension of index x in terms of the cost parame- 
ter X, such that dim(x) = ax^- If legCosts is not 
specified, it is assumed that each index (whether 
summed or unsummed) has dimension %. Note that 
b may take the value zero, permitting a fixed cost 
to be specified for some indices. 

Regardless of the value of costType, if legCosts is spec- 
ified, each index label must appear in the first column 
precisely once. Note that tensors of total dimension 1 
are not supported. 

keepTriv: Setting this value to true suppresses the 
stripping of trivial indices. Setting keepTriv to true is 
not recommended as it may result in poorer performance 
of the Netcon algorithm, and the return of sub-optimal 
contraction sequences for some disjoint tensor networks. 
The only advantage to enabling keepTriv is that it may 
occasionally replace outer products with contraction over 
one or more trivial indices. There is no cost benefit to 
this replacement, but the resulting sequence may look a 
little neater. Default value: false. 

On completion, net con () returns an optimal contrac- 
tion sequence and associated cost. These are specified as 
follows: 

sequence: Sequence over which the indices of the ten- 
sor network should be summed in order to contract the 
network for minimum cost. For the interpretation of this 
sequence, see Appendix IB] 

cost: Specifies the total number of multiplication op- 
erations associated with optimal contraction of the tensor 
network, for example according to the sequence returned 
in sequence. For costType — 2 this value is a number. 
For costType — 1 the cost takes the form of a polynomial 



m X, 



i=0 



(14) 



and this is returned as a 1 x Xmax array whose entries 
cost(i) correspond to the coefficients a^-i. 

Note that when a tensor network involves one or more 
traces, these may always be evaluated before network 
contraction begins. Evaluating a trace involves only ad- 
dition operations, not multiplication, and thus these op- 
erations are relatively cheap and are ignored when com- 
puting the value of cost (though the presence of costs 
associated with tracing over indices will be noted in the 
text output if verbosity > 0). 



C. Sample invocation and output 

As a simple example, consider the tensor network given 
in Fig.[TJ Allowing costType and legCosts to take their 
default values, corresponding to each index having di- 
mension X, netcon may be invoked with verbosity level 1 
by the command 

[sequence cost] =netcon({[-l 1 2 3] , [2 4 5 6] , 

[1 5 7 -3] , [3 8 4 9] , 
[6 9 7 10] , [-2 8 11 12] , 
[10 11 12 -4]},1); 

returning the output 

Looking for solutions of cost D(X~1) 

Looking for solutions of cost D(X~6) 

Looking for solutions of cost D(X~7) 

Looking for solutions of cost D(X~8) 



Best sequence: 
Cost: 



94657213 11 12 8 10 
2X~8 + 2X"7 + 2X~6 + OX'S 
+ OX'4 + 0X~3 + OX'2 + 0X~1 
+ OX'O 



indicating that the cost of contracting this network 
scales as O(x^), and that for any given value of x 
the actual cost of performing this contraction will be 
2x^-(-2x^-l-2x^. The sequence and cost are also returned 
in the variables sequence and cost: 

sequence = [9 4 6 5 7 2 1 3 11 12 8 10] 
cost = [0 0000022 2]. 



IV. PERFORMANCE 

The aim of the Netcon algorithm is to present a prac- 
tical, automated means of evaluating the optimal con- 
traction for a tensor network, as an alternative to human 
evaluation. In order to present a compelling argument 
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for this algorithm, it is necessary to show that an imple- 
mentation of the algorithm is capable of delivering results 
within a timeframe comparable to or better than that of 
human evaluation. It can also be argued that Netcon 
would be attractive even if slightly slower than the hu- 
man, as it would free the human from the necessity of 
performing what is a painstaking and frequently tedious 
process while also eliminating the risk of human error. 

In practice the reference implementation netcon () 
supplied with this paper has proven substantially faster 
than human evaluation for all tensor networks examined 
in this paper, with the additional benefit of providing 
absolute certainty that no more efficient contraction se- 
quence remains to be found. (For more complicated ten- 
sor networks, performing an exhaustive search is not an 
option for a human, and the problem is often instead 
reduced to that of finding the contraction sequence of 
lowest cost possible within a reasonable timeframe.) Per- 
formance data are provided in Table IVl 

Of particular interest is the result for the 4:1 2D MERA 
described in Refs. m and |43. Variational optimisation 



TABLE V. Performance comparison between netconO and 
a human operator determining optimal contraction sequences 
for a selection of tensor networks. For the TEBD algorithm, 
the time quoted is that taken to find the optimal contrac- 
tion sequence for application of an imaginary time evolution 
update [see Figs. 3(i)-(ii) of Ref. [S^l, and has physical dimen- 
sion 2 and bond dimension D. For TTNs it is assumed that 
the Hamiltonian is made up of nearest-neighbour terms on a 
square lattice, and the time is for finding an optimal contrac- 
tion sequence to perform energy minimisation on one tensor 
in the tree. For each MERA algorithm listed, the times given 
are for the network in that algorithm taking the longest time 
to evaluate using netconO. The times ihuman correspond to 
the time for the author to independently identify a contrac- 
tion sequence having the same cost to leading order as that 
given by netconO . All calculations were performed using the 
reference implementation of netconO with costType = 2 and 
with the initial value of xiCap set to 1, running in Matlab 
R2011a on a Macbook Pro with a 2.2GHz Intel Core 17 pro- 
cessor and 8Gb of 1333MHz RAM. The C-|-|~ portion of the 
software was compiled using Apple XCode 4.5.2. 



Tensor network 


Number 
of tensors 


Optimal cost 


t^human 


t^nctcon 


TEBD 


6 


op^) 


49s 


0.04s 


3:1 ID TTN 


5 


o(/) 


10s 


0.04s 


9:1 2D TTN 


8 


o(x^') 


35s 


0.06s 


3:1 ID MERA 


7 


O(X^) 


61s 


0.06s 


2:1 ID MERA 


11 


o(x') 


19m 


0.12s 


9:1 2D MERA 


20 


o(x^«) 


38m 


21s 


4:1 2D MERA 


27 


o(x'') 


a 


4 days 



Human analysis of this network failed to yield a contraction 
sequence cheaper than 0{x^*). See Sec. lIVI of main text for 
discussion. 



of this Ansatz involves the contraction of eleven dis- 
tinct tensor networks (and a further forty-nine which 
are equivalent to reflections or rotations of these eleven) , 
and represents a particularly challenging task for Net- 
con given both the large number of tensors involved (27 
in each diagram) and the large number of interconnec- 
tions between the tensors. Determining optimal contrac- 
tion sequences for each of these networks using netconO 
takes times ranging from one hundred and thirty-one 
minutes to four days, although the eleven different di- 
agrams can all be attacked simultaneously by running 
multiple instances of Matlab across several machines. 
This might seem to be an inconveniently lengthy com- 
putation simply to determine an optimal contraction se- 
quence, but the advantage of running Netcon becomes 
clear when it is noticed that the contraction sequences re- 
turned by netconO have total costs of 0(x^®), whereas 
the cost of optimising this Ansatz is reported in Ref. 
to be of 0(x^^)- Assuming that the value reported in 
Ref. is not a typo, the contraction sequences located 
by the Netcon algorithm are therefore significantly more 
efficient than any previously identified by human evalu- 
ation. 

For simpler networks, it is clear from Table |V] that 
netconO offers a substantial advantage in both speed 
and convenience when it comes to determining optimal 
contraction sequences. For more complex networks, the 
challenge presented both to netconO and to human op- 
erators increases, and the ability of netconO to perform 
an exhaustive search within a practical timeframe con- 
tinues to make it the approach of choice. 



V. DISCUSSION 

This paper has presented a detailed description of the 
Netcon algorithm for determining an optimal sequence 
for the contraction of a tensor network, and the associ- 
ated cost, along with a reference implementation in Mat- 
lab called netconO. It has been shown that this refer- 
ence implementation has a performance which is compet- 
itive in real-world applications, evaluating optimal con- 
traction sequences substantially faster than a human re- 
searcher. By automating this time-consuming and error- 
prone task, the Netcon algorithm is capable of substan- 
tially facilitating the development and implementation of 
more advanced tensor network algorithms and Ansatze. 

Authors who use Netcon, netconO, or any derivative 
works, are requested to cite this paper in resulting pub- 
lications. An example citation might read "Our results 
were computed using an implementation of our varia- 
tional Ansatz in Fortran95, with contraction sequences 
determined by the Netcon algorithm [1]" where "[1]" is 
a citation of this paper. 

The author thanks the Ontario Ministry of Research 
and Innovation ERA for financial support. 
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Appendix A: Preferential nature of pairwise 
contractions 



Appendix B: Interpretation of contraction sequences 
specified as a list of indices 



In Sec. lII ATJ it was stated that an optimal contraction 
sequence for a tensor network may always be realised as 
a series of pairwise contractions. To see this, consider the 
contraction of three tensors, A, B, and C, to yield a single 
tensor D. Let ^ab denote the product of the dimensions 
of all indices on tensor A which connect to tensor B, and 
similarly for ^ac and ^bc- Let S,a denote the product of 
all indices on tensor A which do not connect to either B 
or C, and similarly for ^f, and ^c- The dimension of a set 
containing no indices, is always one. Contracting these 
three tensors as a single process involves a cost of 



For example, for the contraction 



(Al) 
(A2) 



we have 



S,a = dim a, 

?6 = 1, 

^c = dime, 

^ab = dim/3, 

^ac = dim 7, 

and ^bc = dim 5. 



(A3) 



For each element of Df it is necessary to sum over 
£.ab£.ac£.bc different contributions (corresponding to the 
enumeration of indices /?, 7, and S), each involving two 
multiplications, and there are then ^a^b^c entries in Df, 
for the total number of multiplication operations given 
in Eq. (JA1[) . In contrast, pairwise contraction may be 
achieved by any of the sequences ((AB)C), ((AC)B), or 
((BC)A), where ((XY)Z) means "contract tensor X with 
tensor Y, then contract the result with tensor Z". The 
sequence ((AB)C) is readily seen to attract a total cost 
of 



<,at,b^ab<,act,bc i saS&ScSacsi 



be 



(A4) 



multiplication operations, with those for ((AC)B) and 
((BC)A) being achieved by the relevant label permuta- 
tions. Since all parameters in Eqs. (|A1[) and (jA4p take a 
value greater than or equal to one, the value of Eq. (IA4[) 
is always less than or equal to that of Eq. (jAip . The 
argument extends directly to contraction of four or more 
tensors, with the cost of sequential pairwise contraction 
continuing to always be less than or equal to that of more 
complicated contractions, and consequently for any ten- 
sor network it is always possible to identify a minimum- 
cost contraction sequence in which all contractions pro- 
ceed in a pairwise fashion. Note that no assumption has 
been made about the values of the index dimensions, and 
thus this result holds even for contraction sequences in- 
volving outer products, which may be equated with con- 
traction over indices of dimension 1. 



The Netcon algorithm described in this paper takes as 
its input a description of a tensor network where each 
summed index is associated with a positive integer la- 
bel and each open index is associated with a negative 
integer label. As its output the algorithm returns an op- 
timal contraction sequence for the specified tensor net- 
work, specified as a list of positive integer labels possibly 
interspersed with zeros, and the cost of performing this 
contraction (corresponding to the number of multiplica- 
tion operations required). Interpretation of tensor con- 
traction sequences specified in this form is described in 
the paper above, but is summarised in this Appendix for 
convenience. 

Starting with a list of tensors in the network to be 
contracted, and beginning with the first index of the se- 
quence, contraction of a tensor network proceeds as fol- 
lows: 

1; If the sequence list is empty, stop. Contraction of the 

tensor network is complete. 
2: Read the first entry, ii. 
3: If h = 0: 

3a: Read a further x —I entries, for a total of x en- 
tries, denoted ii, . . . ,ix, where x is the largest 
possible value such that all x entries are zero. 

3b: Let n be the number of tensors currently in the 
list of tensors. If a; == n — 1: 

3bl: Using the outer product algorithm given 
below, perform an outer product of all n 
remaining tensors. Denote the result of 
this outer product X. Delete all n tensors 
from the list. Add tensor X to the list. 

3b2: Delete entries ii,...,ix from the se- 
quence. 

3b3: Go to step 1. 
3c: Otherwise (i.e. x ^ n — 1): 

3cl: Delete entries ii,...,ix from the se- 
quence. 

3c2: Read the next y entries from the sequence 
(denoted ji,...,jy), and hst all tensors 
on which indices ji, . . . ,jy appear, where 
y is the largest possible value such that 
all entries ji , . . . , jy are nonzero and the 
number of tensors in the list is precisely 
x + 2. 

3c3: Let these tensors be referred to as 
Ai,...,A^+2. 

3c4: Let Bi and B2 denote the tensors on 
which index ji appears. Identify the 
smallest value of z such that index jz ap- 
pears on a tensor which is neither Bi nor 
B2. Index jz also appears on either Bi or 
B2. Let X denote this tensor (either Bi or 
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B2). Note that tensor X will be a member 
of the list Ai, . . . , Pi.x+2- 

3c5: Using the outer product algorithm given 
below, perform an outer product of all 
tensors Ai , . . . , Kx+2 except for tensor X. 
Let the result of this outer product be de- 
noted Y. 

3c6: Indices ji,...,jy all appear on both X 
and Y. Evaluate the product of ten- 
sors X and Y, denoted (XY), summing 
over all possible configurations of the in- 
dices ji, ... ,jj;. 

3c7: Delete tensors Ai, . . . , A2._|_2 from the list 
of tensors. Add tensor (XY) to the list of 
tensors. 

3c8: Delete indices ji,...,jy from the se- 
quence. 

3c9: Return to step 1. 



Otherwise (i.e. ii 7^ 0), identify which tensors index ii 

appears on. 
If index ii appears on only one tensor, it represents a 

trace. Trace over index zi, delete this index from 

the sequence, and return to step 1. 



6: Otherwise: Index ii appears on two tensors, A and 
B. Read a further x — I indices, for a total of x 
indices, denoted ii, . . . ,ix, where x is the largest 
possible value such that all indices ii, . . . ^i^ appear 
on both tensor A and tensor B. 

7: Evaluate the product of tensors A and B, denoted 
(AB), summing over all possible configurations of 
the indices ii, . . . , i^- 

8: Delete tensors A and B from the list of tensors. Add 
tensor (AB) to the list of tensors. 

9: Delete indices ii,. . . ,1^ from the sequence. 

10: Return to step 1. 

When called upon to perform an outer product of m 
tensors, this should be done by applying the following 
simple algorithm: 

1: Define the total dimension of a tensor as the product 
of the dimensions of its indices. 

2: List the m participating tensors. 

3: Identify the two tensors having the smallest total di- 
mensions. 

4: Remove those two tensors from the list. 

5: Add their outer product to the list. 

6: Repeat steps 3-5 until only one tensor remains. 
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